function [ m_ccc,m_css,m_cww,m_cee,m_cnn,m_cne,m_cnw,m_cse, m_csw ] = rrb_index( Nx , Ny , cc , cs, cw)
%function [ m_ccc,m_css,m_cww,m_cee,m_cnn] = rrb_index( Nx , Ny , cc , cs, cw)

%Construction of RRB Preconditioner
%   Detailed explanation goes here

m_ccc =cc;
m_css =cs;
m_cww =cw;

m_cnn = zeros(Nx+3,Ny+3);
m_cee = zeros(Nx+3,Ny+3);
for jj= 1:Ny+2
    for ii=1:Nx+2
        m_cnn(jj,ii) = -cs (jj+1,ii);
        m_cee (jj,ii) = -cw (jj,ii+1);
    end
end
%generatespy_index(Nx,Ny,m_ccc,m_css,m_cww,m_cee,m_cnn);

m_cne = zeros(Nx+3,Ny+3);
m_cnw = zeros(Nx+3,Ny+3);
m_cse = zeros(Nx+3,Ny+3);
m_csw = zeros(Nx+3,Ny+3);

m_Nx1 = Nx;
m_Nx2 = Ny;

st = 1;
tw = 2;
niv = 0;
while tw < m_Nx2+3
if (niv > 0)
    fac = 1.0 + (st) * 10.0 / ((m_Nx1 + 2) * (m_Nx2 + 2));
    %fac=1.0;
    for jj = st+1:tw:m_Nx2+3
        for ii = 1:tw:m_Nx1+3
            m_ccc(jj,ii) = m_ccc(jj,ii) * fac +m_cne(jj,ii) + m_cnw(jj,ii) +m_cse(jj,ii) + m_csw(jj,ii);
            m_cne(jj,ii)=0;
            m_cnw(jj,ii)=0;
            m_cse(jj,ii)=0;
            m_csw(jj,ii)=0;
        end
    end
    
    for jj =1:tw:m_Nx2+3
        for ii = st+1:tw:m_Nx1+3
            m_ccc(jj,ii) = m_ccc(jj,ii) * fac +m_cne(jj,ii) + m_cnw(jj,ii) +m_cse(jj,ii) + m_csw(jj,ii);
            m_cne(jj,ii)=0;
            m_cnw(jj,ii)=0;
            m_cse(jj,ii)=0;
            m_csw(jj,ii)=0;
        end
    end
end

niv=niv+1;
% Gaussian Elimination: This is done for all R1 and R2 Nodes.

for jj = st+1:tw:m_Nx2+3
    for ii = st+1:tw:m_Nx1+3
        iw = ii - st;
        ie = ii + st;
        js = jj - st;
        jn = jj + st;
        if (jn <= m_Nx2 + 3)
            coeff = m_cnn(jj,ii)/ m_ccc(jn,ii);
            m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jn,ii)*coeff;
            m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jn,ii)* coeff;
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jn,ii) *coeff;
            m_cnn(jj,ii) = -m_cnn(jn,ii) *coeff;
        end
        if (js > 1)
            
            m_cse(jj,ii) = m_cse(jj,ii) - m_css(jj,ii) * m_cee(js,ii) / m_ccc(js,ii);
            m_csw(jj,ii) = m_csw(jj,ii) - m_css(jj,ii) * m_cww(js,ii) / m_ccc(js,ii);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jj,ii) * m_cnn(js,ii) / m_ccc(js,ii);
            m_css(jj,ii) = -m_css(jj,ii) * m_css(js,ii) / m_ccc(js,ii);
        end
        if (ie <= m_Nx1 + 3)
            
            m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jj,ii) * m_cnn(jj,ie) / m_ccc(jj,ie);
            m_cse(jj,ii) = m_cse(jj,ii) - m_cee(jj,ii) * m_css(jj,ie) / m_ccc(jj,ie);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cee(jj,ii) * m_cww(jj,ie) / m_ccc(jj,ie);
            m_cee(jj,ii) = -m_cee(jj,ii) * m_cee(jj,ie) / m_ccc(jj,ie);
        end
        if (iw > 1)
            
            m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jj,ii) * m_cnn(jj,iw) / m_ccc(jj,iw);
            m_csw(jj,ii) = m_csw(jj,ii) - m_cww(jj,ii) * m_css(jj,iw) / m_ccc(jj,iw);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cww(jj,ii) * m_cee(jj,iw) / m_ccc(jj,iw);
            m_cww(jj,ii) = -m_cww(jj,ii) * m_cww(jj,iw) / m_ccc(jj,iw);
        end
    end
end

for jj = 1:tw:m_Nx2+3
    for ii = 1:tw:m_Nx1+3
        iw = ii - st;
        ie = ii + st;
        js = jj - st;
        jn = jj + st;
        if (jn <= m_Nx2 + 3)
            m_cne(jj,ii) = m_cne(jj,ii) - m_cnn(jj,ii) * m_cee(jn,ii) / m_ccc(jn,ii);
            m_cnw(jj,ii) = m_cnw(jj,ii) - m_cnn(jj,ii) * m_cww(jn,ii) / m_ccc(jn,ii);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cnn(jj,ii) * m_css(jn,ii) / m_ccc(jn,ii);
            m_cnn(jj,ii) = -m_cnn(jj,ii) * m_cnn(jn,ii) / m_ccc(jn,ii);
        end
        if (js > 1)
            
            m_cse(jj,ii) = m_cse(jj,ii) - m_css(jj,ii) * m_cee(js,ii) / m_ccc(js,ii);
            m_csw(jj,ii) = m_csw(jj,ii) - m_css(jj,ii) * m_cww(js,ii) / m_ccc(js,ii);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jj,ii) * m_cnn(js,ii) / m_ccc(js,ii);
            m_css(jj,ii) = -m_css(jj,ii) * m_css(js,ii) / m_ccc(js,ii);
        end
        if (ie <= m_Nx1 + 3)
            
            m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jj,ii) * m_cnn(jj,ie) / m_ccc(jj,ie);
            m_cse(jj,ii) = m_cse(jj,ii) - m_cee(jj,ii) * m_css(jj,ie) / m_ccc(jj,ie);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cee(jj,ii) * m_cww(jj,ie) / m_ccc(jj,ie);
            m_cee(jj,ii) = -m_cee(jj,ii) * m_cee(jj,ie) / m_ccc(jj,ie);
        end
        if (iw > 1)
            
            m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jj,ii) * m_cnn(jj,iw) / m_ccc(jj,iw);
            m_csw(jj,ii) = m_csw(jj,ii) - m_cww(jj,ii) * m_css(jj,iw) / m_ccc(jj,iw);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cww(jj,ii) * m_cee(jj,iw) / m_ccc(jj,iw);
            m_cww(jj,ii) = -m_cww(jj,ii) * m_cww(jj,iw) / m_ccc(jj,iw);
        end
    end
end

for jj = st+1:tw:m_Nx2+3
    for ii = st+1:tw:m_Nx1+3
        m_ccc(jj,ii) = m_ccc(jj,ii) +m_cnn(jj,ii) +m_cee(jj,ii) +m_css(jj,ii) +m_cww(jj,ii);
        
        %                 // Om geheugen te besparen gaan we cnw(jj,ii) opslaan in cnn(jj,ii);
        %                 //                                 cne(jj,ii) opslaan in cee(jj,ii);
        %                 //                                 csw(jj,ii) opslaan in cww(jj,ii);
%         %                 //                                 cse(jj,ii) opslaan in css(jj,ii);
%         m_cnn(jj,ii) = m_cnw(jj,ii);
%         m_cee(jj,ii) = m_cne(jj,ii);
%         m_cww(jj,ii) = m_csw(jj,ii);
%         m_css(jj,ii) = m_cse(jj,ii);
        m_cnn(jj,ii) = 0;
        m_cee(jj,ii) = 0;
        m_cww(jj,ii) = 0;
        m_css(jj,ii) = 0;
        
    end
end
niv=niv+1;

%
%             // Uitgaande van het stencil
%         //                          cnn
%         //                     cnw      cne
%         //                 cww      ccc      cee
%         //                     csw      cse
%         //                          css
%         // elimineren we cnw en cne met de scheve 5-punts formule. Dit geeft
%         //                 cnw      cnn      cne
%         //                 cww      ccc      cee
%         //                 csw      css      cse
%         // We zijn nu weer terug de eerste 9-punts formule, maar op een
%         // grover rooster.
for jj = 1:tw:m_Nx2+3
    for ii = 1:tw:m_Nx1+3
        iw = ii - st;
        ie = ii + st;
        js = jj - st;
        jn = jj + st;
        if (jn <= m_Nx2 + 3 && ie <= m_Nx1 + 3)
            
            m_cnn(jj,ii) = m_cnn(jj,ii) - m_cne(jj,ii) * m_cnw(jn,ie) / m_ccc(jn,ie);
            m_cee(jj,ii) = m_cee(jj,ii) - m_cne(jj,ii) * m_cse(jn,ie) / m_ccc(jn,ie);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cne(jj,ii) * m_csw(jn,ie) / m_ccc(jn,ie);
            m_cne(jj,ii) = -m_cne(jj,ii) * m_cne(jn,ie) / m_ccc(jn,ie);
        end
        if (iw > 1 && jn <= m_Nx2 + 3)
            
            m_cnn(jj,ii) = m_cnn(jj,ii) - m_cnw(jj,ii) * m_cne(jn,iw) / m_ccc(jn,iw);
            m_cww(jj,ii) = m_cww(jj,ii) - m_cnw(jj,ii) * m_csw(jn,iw) / m_ccc(jn,iw);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cnw(jj,ii) * m_cse(jn,iw) / m_ccc(jn,iw);
            m_cnw(jj,ii) = -m_cnw(jj,ii) * m_cnw(jn,iw) / m_ccc(jn,iw);
        end
        if (js > 1 && iw > 1)
            
            m_css(jj,ii) = m_css(jj,ii) - m_csw(jj,ii) * m_cse(js,iw) / m_ccc(js,iw);
            m_cww(jj,ii) = m_cww(jj,ii) - m_csw(jj,ii) * m_cnw(js,iw) / m_ccc(js,iw);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_csw(jj,ii) * m_cne(js,iw) / m_ccc(js,iw);
            m_csw(jj,ii) = -m_csw(jj,ii) * m_csw(js,iw) / m_ccc(js,iw);
        end
        if (js > 1 && ie <= m_Nx1 + 3)
            
            m_css(jj,ii) = m_css(jj,ii) - m_cse(jj,ii) * m_csw(js,ie) / m_ccc(js,ie);
            m_cee(jj,ii) = m_cee(jj,ii) - m_cse(jj,ii) * m_cne(js,ie) / m_ccc(js,ie);
            m_ccc(jj,ii) = m_ccc(jj,ii) - m_cse(jj,ii) * m_cnw(js,ie) / m_ccc(js,ie);
            m_cse(jj,ii) = -m_cse(jj,ii) * m_cse(js,ie) / m_ccc(js,ie);
        end
    end
end
st = tw;
tw = tw*2;

% for jj = 1:1:m_Nx2+3
%     for ii = 1:1:m_Nx1+3
%         
%         m_ccc(jj,ii) = 1.0 / m_ccc(jj,ii);
%         m_cnn(jj,ii) = m_cnn(jj,ii) * m_ccc(jj,ii);
%         m_css(jj,ii) = m_css(jj,ii) * m_ccc(jj,ii);
%         m_cee(jj,ii)= m_cee(jj,ii) * m_ccc(jj,ii);
%         m_cww(jj,ii) = m_cww(jj,ii) * m_ccc(jj,ii);
%         m_cnw(jj,ii) = m_cnw(jj,ii) * m_ccc(jj,ii);
%         m_cne(jj,ii) = m_cne(jj,ii) * m_ccc(jj,ii);
%         m_csw(jj,ii) = m_csw(jj,ii) * m_ccc(jj,ii);
%         m_cse(jj,ii) = m_cse(jj,ii) * m_ccc(jj,ii);
%        
%         
%     end
% end

end


end

